	subroutine simple_diir2d( nfx,nfy, ngx, ngy, f, g, h)
	double precision f(0:(nfx-1),0:(nfy-1))
	double precision g(0:(ngx-1),0:(ngy-1))
	double precision h(0:(nfx-1),0:(nfy-1))
	integer nfx, nfy, ngx, ngy

	integer i, j, ih, jh, i1, i2, j1, j2

	do j = 0, nfy-1
	    do i = 0, nfx-1
		h(i,j) = f(i,j)
	    end do
	end do

	do jh = 0, nfy-1
	    do ih = 0, nfx-1
		j1 = max( 0, jh - (ngy-1))
		j2 = jh - 1
		i1 = max( 0, ih - (ngx-1))
		i2 = ih - 1
		do j = j1, j2
		    do i = i1, ih
			h(ih,jh) = h(ih,jh) - h(i,j) * g(ih-i, jh-j)
		    end do
		end do
		    do i = i1, i2
			h(ih,jh) = h(ih,jh) - h(i,jh) * g(ih-i, 0)
		    end do
		    h(ih,jh) = h(ih,jh) / g(0,0)
	    end do
	end do
	return
	end
